/*************************** model.cpp **********************************
* Author:        Gena Tang
* Date created:  2012-09-12
*****************************************************************************/
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <stdio.h>
#include <map>
#include <sstream>
#include <string.h>
#include <algorithm>
#include <math.h>

#include "model.h"     // class definition

Gene::Gene(string line)
{
        size_t index = 0;
        while (true) {    
     index = line.find(">", index);
     if (index == string::npos) break;
     line.replace(index, 1, " ");
     ++index;
        }
        index = 0;
        while (true) {    
     index = line.find("<", index);
     if (index == string::npos) break;
     line.replace(index, 1, " ");
     ++index;
        }

        istringstream str(line);
        str >> id >> strand >> gene_left >> gene_right >> gene_length >> gene_class;
//      cout <<"initialize: geneleft="<< gene_left << " gene_right="<< gene_right << " gene_class="<< gene_class<<endl;
//if(id == 1831){
//              cout << endl;
//      }
};

void Gene::seqToInt(string & str){
        replace( str.begin(), str.end(), 'A', '0');
        replace( str.begin(), str.end(), 'C', '1');
        replace( str.begin(), str.end(), 'G', '2');
        replace( str.begin(), str.end(), 'T', '3');
};

void Gene::set_gene_seq(string & seq)
{
        seq.erase(remove(seq.begin(),seq.end(),'\n'),seq.end());
        seq.erase(remove(seq.begin(),seq.end(),'\r'),seq.end());
        if(strand == '+'){

                gene_seq.assign( seq.substr(gene_left - 1, gene_length) );
        //      cout <<" seq="<< gene_seq<<endl;
        //      cout <<" left "<< gene_left<<endl;
                //cout <<"set gene_seq left=" << gene_left << " right="<< gene_right << " seq="<< gene_seq<<endl;
                seqToInt(gene_seq);
        }
        else
        {
                gene_seq.assign( revcomp( seq.substr(gene_left - 1, gene_length)  ) );
                //cout << "set gene_seq left="<<gene_left << " right="<< gene_right << " seq="<< gene_seq<<endl;
                seqToInt(gene_seq);
        }
};

void Gene::set_upstream_seq(string & seq, int l)
{
        if(strand == '+'){
                if(gene_left > l + 1)
                        upstream_seq.assign( seq.substr(gene_left - 1 - l, l) );
                else
                        upstream_seq.assign( seq.substr(0, l) );
                seqToInt(upstream_seq);
        }
        else
        {
                if(gene_right < seq.length() - l )
                        upstream_seq.assign( revcomp( seq.substr(gene_right, l)  ) );
                else
                        upstream_seq.assign( revcomp( seq.substr(gene_right, seq.length()-gene_right)  ) );
        }
};

double Gene::get_codon_score(double * P_Q){
        double score = 0;
        for( int i = 3; i < gene_seq.size() - 4; i += 3){
                int index = 0;
                index += 16*atoi( gene_seq.substr(i,1).c_str() );
                index += 4*atoi( gene_seq.substr(i+1,1).c_str() );
                index += atoi( gene_seq.substr(i+2,1).c_str() );
                if(index == 49)
                        index = 48;
                else if(index > 56)
                        index -= 3;
                else if(index > 49)
                        index -= 2;
                score += P_Q[index];
        }
        return score;
};

double Gene::get_start_score(double * P_Q){

        string x = gene_seq.substr(0,1);
        //cout << x << " " << gene_seq <<endl;
        if(x == "0")
                return P_Q[0];
        if(x == "2")
                return P_Q[1];
        if(x == "3")
                return P_Q[2];

        return 0.0;

};
double Gene::get_length_score(vector<double>& P_Q){

        int size = gene_seq.size()/3-1; //exclude stop codon
        if( size <= P_Q.size() )
                return P_Q[size-1];

        return 0.0;

};



void Sequence::read(string name)
{
        ifstream file;
        file.open( name.c_str() );
        bool SEQNAME = false;
        if ( !file.is_open() )
        {
                cout << "error on open file: " << name << endl;
                cin.get();
        }
        string line;
        double value;

        while( !file.eof() )
        {
                getline( file, line );

                if ( line.empty() )
                        continue;
                if ( line.find_first_of('>') != string::npos ){
                        name.assign( line.substr( 1, line.size()-1 ) );
                        SEQNAME=true;
                        continue;
                }
                if(SEQNAME){
                        seq.append( line );
                }
        }
};

Sequence::Sequence(){
} 



ModelFile::ModelFile(){
};
/***********************************************************************
read model file function, default order is 2
***********************************************************************/
void ModelFile::read(string name)
{
        // test if order is 2, if not, exit
        ifstream file;

        file.open( name.c_str() );
        if ( !file.is_open() )
        {
                cout << "error on open file: " << name << endl;
                cin.get();
        }

        


        bool COD1 = false;
        bool COD2 = false;
        bool NON = false;
        bool RBSM = false;
        bool RBSD = false;
        bool UseModelDur = true;
		bool STARTM = false;

        string line;
        double value;

        while( !file.eof() )
        {
                getline( file, line );

                if ( line.empty() )
                        continue;

                if ( line.at(0) == '#' )
                        continue;


                if (  ! line.compare( 0, 4, "ORDM" ) ){
                        order = atof(line.substr(5).c_str());
                        if( order != 2){
                                cout << "error on model file: model order is not 2" << endl;
                                cin.get();
                        }
						msize =  4 << (2*order);
						P_1.reserve( msize );
						H_1.reserve( msize );
						Q.reserve( msize );
						P_2.reserve( msize );
						H_2.reserve( msize );
						P_3.reserve( msize );
						H_3.reserve( msize );
						P_1_rel.reserve( msize );
						H_1_rel.reserve( msize );
						Q_rel.reserve( msize );
						P_2_rel.reserve( msize );
						H_2_rel.reserve( msize );
						P_3_rel.reserve( msize );
						H_3_rel.reserve( msize );
						P.reserve( msize );
						for(int i = 0; i < msize; i ++){
							P_1_rel.push_back(0.0);
							H_1_rel.push_back(0.0);
							Q_rel.push_back(0.0);
							P_2_rel.push_back(0.0);
							H_2_rel.push_back(0.0);
							P_3_rel.push_back(0.0);
							H_3_rel.push_back(0.0);
							P.push_back(0.0);
						}
                }
				if (  ! line.compare( 0, 4, "NDEC" ) ){
                        noncodDecay = atof(line.substr(5).c_str());
                        
                }
				if (  ! line.compare( 0, 4, "CDEC" ) ){
                        codDecay = atof(line.substr(5).c_str());
                       
                }
				if (  ! line.compare( 0, 4, "CDCD" ) ){
                        codToCod = atof(line.substr(5).c_str());
                      
                }
                if (  ! line.compare( 0, 4, "RBSS" ) ){
                        motif_l = atof(line.substr(5).c_str());
                        cout << "motif length is" << motif_l << endl;
                }
                if (  ! line.compare( 0, 4, "RBSL" ) ){
                        spacer_l = atof(line.substr(5).c_str());
                        cout << "spacer length is " << spacer_l << endl;
                }
                if (  ! line.compare( 1, 2, "TG" ) ){
                        P_start.push_back(atof(line.substr(5).c_str()));
						if( atof(line.substr(5).c_str()) != 0 ){
							start_codons.push_back(line.substr(0,3));
						}
                }
				if (  ! line.compare( 0, 3, "TAG" )  ||  ! line.compare( 0, 3, "TGA" ) ||  ! line.compare( 0, 3, "TAA" )){
						if( atof(line.substr(5).c_str()) != 0 ){
							stop_codons.push_back(line.substr(0,3));
						}
                }
				if (  ! line.compare( 0, 4, "MINC" ) ){
                        minc = atof(line.substr(5).c_str());
                        cout << "min coding length is " << minc << endl;
                }
				if (  ! line.compare( 0, 4, "MAXC" ) ){
                        maxc = atof(line.substr(5).c_str());
                        cout << "max coding length is " << maxc << endl;
                }
				if (  ! line.compare( 0, 4, "MAXN" ) ){
                        maxn = atof(line.substr(5).c_str());
                        cout << "mas noncoding length is " << maxn << endl;
                }
                if ( ! isdigit( line.at(0) ))
                {
                        COD1 = false;
                        COD2 = false;
                        NON  = false;
                        RBSM = false;
                        RBSD = false;
						STARTM = false;

                        if (  ! line.compare( 0, 4, "COD1" ) )
                                COD1 = true;
                        if (  ! line.compare( 0, 4, "COD2" ) )
                                COD2 = true;
                        if (  ! line.compare( 0, 4, "NONC" ) )
                                NON  = true;
                        if (  ! line.compare( 0, 4, "RBSM" ) )
                                RBSM  = true;
                        if (  ! line.compare( 0, 4, "RBSD" ) )
                                RBSD  = true;
						if (  ! line.compare( 0, 6, "STARTM" ) )
                                STARTM  = true;
                        continue;
                }

                value =  atof( line.c_str() );

                if ( value < 0 )
                {
                        cout << "error: negative probability value found in file " << name << " : " << line <<  endl;
                        cin.get();
                }

                if ( NON )
                        Q.push_back( value );
                if ( COD1 ){
					double frame1, frame2, frame3;
					istringstream str(line);
                    str >> frame1 >> frame2 >> frame3;
                    P_1.push_back( frame1 );
					P_2.push_back( frame2 );
					P_3.push_back( frame3 );
				}
                if ( COD2 ){
                    double frame1, frame2, frame3;
					istringstream str(line);
                    str >> frame1 >> frame2 >> frame3;
                    H_1.push_back( frame1 );
					H_2.push_back( frame2 );
					H_3.push_back( frame3 );
				}
                if ( RBSD )
                        spacer.push_back( value );
                if ( RBSM )
                {
                        vector<double> vec; // composition at on position of a motif
                        double Ant,Cnt,Gnt,Tnt;
                        istringstream str(line);
                        str >> Ant >> Cnt >> Gnt >> Tnt;
                        vec.push_back(Ant);
                        vec.push_back(Cnt);
                        vec.push_back(Gnt);
                        vec.push_back(Tnt);
                        RBS.insert(RBS.begin(), vec);

                }
				if ( STARTM )
                {
                        vector<double> vec; // composition at on position of a motif
                        double Ant,Cnt,Gnt,Tnt;
                        istringstream str(line);
                        str >> Ant >> Cnt >> Gnt >> Tnt;
                        vec.push_back(Ant);
                        vec.push_back(Cnt);
                        vec.push_back(Gnt);
                        vec.push_back(Tnt);
                        StartM.push_back(vec);

                }
        }

        file.close();

        if ( P_1.size() != 64 && Q.size() != 64 )
        {
                cout << "error: P or Q arrays are of unexpected size. only work for second order right now" << endl;
                cin.get();
        }


        int ATG = ((A<<4)+(T<<2)+G);
        int CTG = ((C<<4)+(T<<2)+G);
        int GTG = ((G<<4)+(T<<2)+G);
        int TTG = ((T<<4)+(T<<2)+G);



        if( UseModelDur )
                dur_non = durationFromModel( Q );

        cout << "NON_DUR " << dur_non << endl;

        getZero( Q, Q_4 );
        normArr( Q_4 );
		normArr(P_start);
		getCodLengthDistr();
		getNonLengthDistr();

        //printVector("Q_4", Q_4 );
	/*	AbsToFirst( P_1, P_1_rel );
		AbsToFirst( P_2, P_2_rel );
		AbsToFirst( P_3, P_3_rel );
		AbsToFirst( H_1, P_1_rel );
		AbsToFirst( H_2, H_2_rel );
		AbsToFirst( H_3, H_3_rel );
		AbsToFirst( Q, Q_rel );*/

		AbsToLast( P_1);
		AbsToLast( P_2 );
		AbsToLast( P_3);
		AbsToLast( H_1 );
		AbsToLast( H_2 );
		AbsToLast( H_3 );
		AbsToLast( Q );
		avg( P, P_1, P_2, P_3);
		
};

void ModelFile::getCodLengthDistr(){
	double temp;
	//i is gene length in nt. codLengthDistr[0] corresponds to minlength
	 for (double i = minc; i <= maxc; ++i )
        {
                if ( (int)i%3 == 0 )
                {
                        temp =  (1/(2*codDecay*codDecay*codDecay)) * i * i * exp( - (double)i / codDecay ) ;
						codLengthDistr.push_back(temp);
                } else
                        codLengthDistr.push_back(0);
        }
	 normArr( codLengthDistr );
	  for (int i = 0; i < codLengthDistr.size(); ++i ){
		  if(codLengthDistr[i] <= 0 ){
			  codLengthDistr[i] = -999;
		  }
		  else {codLengthDistr[i] = log( codLengthDistr[i] );}
	  }

}

void ModelFile::getNonLengthDistr(){
		double temp;
        for (double  i = 1; i <= maxn; ++i ) {
                temp= (1/noncodDecay) * exp( - i / noncodDecay );
				nonLengthDistr.push_back(temp);
        }
		normArr( nonLengthDistr );
		for (int i = 0; i < nonLengthDistr.size(); ++i ){
		  nonLengthDistr[i] = log( nonLengthDistr[i] );
	  }
}

void ModelFile::getZero( vector<double>& arr, vector<double>& outarr )
{

        if(arr.size() != 64)
        {
                cout << "Input is not second order!" << endl;
                cin.get();
        }
        outarr.resize(4,0);
        for(int i = 0; i < 4; ++ i)
        {
                for( int j = 0; j < 4; ++ j)
                {
                        for( int k = 0; k < 4; ++ k)
                        {
                                int pos = ((i<<4)+(j<<2)+k);
                                outarr[i] += arr[pos];
                                outarr[j] += arr[pos];
                                outarr[k] += arr[pos];
                        }
                }
        }
};
// ------------------------------------
double ModelFile::durationFromModel( vector<double>& arr )
{
        int A = 0;
        int C = 1;
        int G = 2;
        int T = 3;

        int TAA = ((T<<4)+(A<<2)+A);
        int TAG = ((T<<4)+(A<<2)+G);
        int TGA = ((T<<4)+(G<<2)+A);

        double x = arr[TAA] + arr[TAG] + arr[TGA];
        x = 1/x;

        return x;
};

// ------------------------------------
void ModelFile::removeStopCodons( vector<double>& arr )
{
        int A = 0;
        int C = 1;
        int G = 2;
        int T = 3;

        int TAA = ((T<<4)+(A<<2)+A);
        int TAG = ((T<<4)+(A<<2)+G);
        int TGA = ((T<<4)+(G<<2)+A);

        // order is important
        arr.erase( arr.begin() + TGA );
        arr.erase( arr.begin() + TAG );
        arr.erase( arr.begin() + TAA );


};

// ------------------------------------
void ModelFile::normArr( vector<double>& arr )
{
        double sum = 0;

        for( int i = 0; i < arr.size(); ++i )
                sum += arr[i];

        if ( !sum )
        {
                cout << "error: sum of probailities is zero" << endl;
                cin.get();
        }

        for( int i = 0; i < arr.size(); ++i )
                arr[i] /= sum;

};

void ModelFile::AbsToFirst( vector<double> & abs, vector<double>& rel )
{
	double temp;

	const int leftA = ( A << (2*order) );
	const int leftC = ( C << (2*order) );
	const int leftG = ( G << (2*order) );
	const int leftT = ( T << (2*order) );
	
	int i = 0;
	while ( i < msize/4 ) {

		temp = abs[leftA] + abs[leftC] + abs[leftG] + abs[leftT];

		if ( temp != 0 ) {
			rel[leftA+i]= abs[leftA+i] / temp;
			rel[leftC+i]= abs[leftC+i] / temp;
			rel[leftG+i]= abs[leftG+i] / temp;
			rel[leftT+i]= abs[leftT+i] / temp;
		} else {
			rel[leftA+i] = 0;
			rel[leftC+i] = 0;
			rel[leftG+i] = 0;
			rel[leftT+i] = 0;
		}

		i++;
	}
};

void ModelFile::AbsToLast( vector<double> & abs )
{
      
        double temp;
		int i = 0;
        while ( i < abs.size() ) {

                temp = abs[A + i] + abs[C + i] + abs[G + i] + abs[T + i];

                if ( temp != 0 ) {
                        abs[A + i] /= temp;
                        abs[C + i] /= temp;
                        abs[G + i] /= temp;
                        abs[T + i] /= temp;
                } else {
                        abs[A + i] = 0;
                        abs[C + i] = 0;
                        abs[G + i] = 0;
                        abs[T + i] = 0;
                }

                i += 4;
        }
};

void printVector(string name, vector<double>& arr)
{
        cout << name << "("<< arr.size() << ")=" << endl;
        for(int i = 0; i < arr.size(); ++i)
        {
                cout << arr[i] << endl;
        }
};


string revcomp(string str){
        reverse( str.begin(), str.end() );
        replace( str.begin(), str.end(), 'A', '3');
        replace( str.begin(), str.end(), 'C', '2');
        replace( str.begin(), str.end(), 'G', '1');
        replace( str.begin(), str.end(), 'T', '0');
        return str;
};

string revcompToNT(string str){
        string s = revcomp(str);
		replace( s.begin(), s.end(), '0', 'A');
        replace( s.begin(), s.end(), '1', 'C');
        replace( s.begin(), s.end(), '2', 'G');
        replace( s.begin(), s.end(), '3', 'T');

        return s;
};

// ------------------------------------
double* VectorToArray( vector<double>& arr )
{
        double* ptr = new double[ arr.size() ];
        if ( !ptr ) 
        {
                cout << "error: out of memory" << endl;
                cin.get();
        }

        for( int i = 0; i < arr.size(); ++i )
                ptr[i] = arr[i];

        return ptr;
};
// ------------------------------------
int* VectorToArray( vector<int>& arr )
{
        int* ptr = new int[ arr.size() ];
        if ( !ptr ) 
        {
                cout << "error: out of memory" << endl;
                cin.get();
        }

        for( int i = 0; i < arr.size(); ++i )
                ptr[i] = arr[i];

        return ptr;
};

// ------------------------------------
double* LogOdd( double* F, double* S, int size )
{
        double* ptr = new double[ size ];
        if ( !ptr ) 
        {
                cout << "error: out of memory" << endl;
                cin.get();
        }

        for( int i = 0; i < size; ++i )
        {
                if ( F[i] <= 0 || S[i] <= 0 )
                {
                        cout << "error: wrong value was found under log: " << F[i] << " "  << S[i] << endl;
                        cin.get();
                }

                ptr[i] = log( F[i]/S[i] );
        }


        return ptr;
};

// ------------------------------------
void LengthLogOdd( vector<double>& logodd, int min, int max, double gamma, double exp )
{
        logodd.resize( max, 0 );
        for( int i = 0; i < max; ++i )
                logodd[i] = 0;

        double P_dur_norm = NormGammConstant( gamma, min, max );
        double Q_dur_norm = NormExpConstant ( exp,   min, max );

        double x;

        for( int i = min; i <= max; ++i )
        {
                x = GetGamma( i, gamma ) / P_dur_norm;
                x /= ( GetExp( i, exp )  / Q_dur_norm );

                logodd[ i - 1 ] = log(x);

        }
};

// ------------------------------------
double GetExp( int i, double betta )
{
        return (1.0/betta) * exp( -i / betta );
};

// ------------------------------------
double GetGamma( int i, double betta )
{
        return (1.0/(2*betta*betta*betta)) * i * i * exp( -i / betta );
};

// ------------------------------------
double NormExpConstant( double betta, int min, int max )
{
        double x = 0;

        for( int i = min; i <= max; ++i )
                x += GetExp( i, betta );


        return x;
};

// ------------------------------------
double NormGammConstant( double betta, int min, int max )
{
        double x = 0;

        for( int i = min; i <= max; ++i )
                x += GetGamma( i, betta );


        return x;
};

void avg(vector<double>& p, vector<double>& p1, vector<double>& p2, vector<double>& p3){
	if(p.size() != p1.size()){
		cout << "error!"<< endl;
		cin.get();
		exit(1);
	}
	for(int i = 0; i < p.size(); ++ i){
		p[i] = 1.0/3 * (p1[i] + p2[i] + p3[i]);
	}
}